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Abstract. 

We examine the conditions under which binary and multiple stars 
may form out of turbulent molecular cloud cores using high resolution 
3-D, adaptive mesh refinement (AMR) hydrodynamics (Truelove et al., 
1997, 1998; Klein, 1999). We argue that previous conclusions on the con- 
ditions for cloud fragmentation have limited applicability, since they did 
not allow for the nonlinear density and velocity perturbations that are 
ubiquitous in molecular cloud cores. Over the past year, we have begun to 
simulate the evolution of marginally stable, turbulent cores. These mod- 
els have radii, masses, density contrasts, turbulent linewidths, and pro- 
jected velocity gradients consistent with observations of low-mass molec- 
ular cloud cores. Our models are evolved in time under self-gravitational 
hydrodynamics with AMR using a barotropic equation of state that mod- 
els the transition from an isothermal to an adiabatic equation of state. 
We examine several properties of the resulting protostellar fragments and 
discuss the qualitative nature of the fragmentation process in realistic 
cloud cores: the transition from single to binary and multiple stars; the 
formation of misaligned binary systems; and the role played by filament 
formation in the formation of stars. 



1. Introduction 



Most stars exist in gravitationally bound binary and low-order multiple systems 
(Bodenheimer et al 2000). Moreover, the observed trend of decreasing fraction 
of binary and multiple systems with age strongly suggests that multiple star 
systems must have formed via fragmentation during the earliest stages of cloud 
collapse, rather than from capture of two or more stars formed individually (i.e. 
Simon et al. 1995). Multiple star formation is central to our understanding 
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of the star formation process, and has enormous significance for the study of 
formation of planetary systems as well. 

While the results of previous authors (notably Hoylc 1953, Hunter 1967, 
Low k Lynden-Bell 1977, Silk 1977, Rees 1977, Larson 1985, and Inutsuka k 
Miyama 1992, among many others) have advanced the state of our knowledge of 
the dynamics underlying the fragmentation process, no theoretical fragmentation 
models to date have been able to convincingly explain the observed properties 
of binaries. Most analytic and semi-analytic work has treated only the earliest 
phases of the collapse, typically via a linearized analysis, and most numerical 
work has treated only the earliest isothermal phase of the collapse, well before 
most of the mass of the cloud core has been accreted onto the fragments. In 
addition one must also be able to follow the collapse long enough to determine 
the number of protostars, their masses, and their orbits. Long before this point, 
the gas ceases to be isothermal. Hence, a nonlinear, nonisothermal analysis is 
absolutely essential. 

Until very recently, the extreme variations in length scale inherent in col- 
lapse have made it difficult to perform accurate calculations of the fragmentation 
process, which is intrinsically three-dimensional in nature. Our recent develop- 
ment of a robust adaptive mesh refinement (AMR) self-gravitational hydrody- 
namics code has resulted in our discovery of a physically motivated constraint 
on the resolution (the Jeans condition) that must be satisfied to suppress arti- 
ficial fragmentation (Truelove et al. 1997). This condition indicates that many 
published solutions were underrcsolved and therefore gave misleading results. 
Since our code is able to adhere to the Jeans criterion, it is ideally suited for the 
treatment of the collapse and fragmentation of molecular cloud cores. 

Over the last year, we have begun a program of research investigating 
the properties of marginally stable, turbulent low-mass molecular cloud cores. 
These models have radii, masses, density contrasts, turbulent linewidths, and 
projected velocity gradients consistent with observations of low-mass molecular 
cloud cores. 

1.1. Key Questions to be Addressed 

Answers to basic questions concerning the fragmentation process remain elusive: 

(1) Why do most stars form in binary systems? 

(2) What determines the efficiency of star formation, i.e., the fraction of a 
core that will end up in protostellar objects? 

(3) What are the properties of protostellar disks formed in turbulent molec- 
ular cloud cores? Are they typically aligned or misaligned relative to one an- 
other? 

(4) What determines the distribution of fragment masses, which is related 
to the initial distribution of stellar masses? Similarly, what determines the initial 
orbital angular momenta and rotational spins of the individual fragments? 

Our goal is to make significant advances in our understanding of all of these 
questions. We shall study the collapse and fragmentation of of turbulent molec- 
ular cloud cores with conditions typical of nearby star-forming regions (Jijina et 
al. 1999), and determine whether our models are consistent with the observed 
frequency of binary star systems. To the best of our knowledge, no previous 
study of binary formation has considered either realistic near-equilibrium clouds 
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or realistic perturbation spectra. This study has the potential of leading to 
a major advance in our understanding of how fragmentation proceeds in real 
molecular cloud cores. 



2. The Key Role of Nonlinear Perturbations : A Numerical Exper- 
iment 

In order to reach more general conclusions on the outcome of the collapse, sev- 
eral groups have attempted to establish a general criterion as to whether an 
individual cloud model will fragment or not. Several classic studies were con- 
ducted based on the results of early axisymmctric and 3D numerical models. 
Notably, Tohline, Durisen, and McCollough (1985), Hachisu & Eriguchi (1985) 
and Miyama (1992) assumed that the collapse of rigidly rotating, uniform den- 
sity spheres resulted in a self-gravitating, rotating equilibrium isothermal disks 
consisting of the entire mass of the cloud. By analyzing these disks for linear 
stability to growing modes using the analysis of Goldreich & Lynden-Bell (1965), 
they arrived at the classic fragmentation criterion for uniform cores that is still 
frequently cited in the literature today. Their criterion is stated in terms of 
the standard two-dimensional parameter space of the cores: a = 5clRo/2GM, 
which is proportional to the ratio of thermal to gravitational energies within the 
core, and /3 = Q,qRq/3GM, which is proportional to the ratio of rotational to 
gravitational binding energies within the core. In terms of a and /?, the classic 
criterion for fragmentation is a/? < 0.12, suggesting that typical slowly rotating 
cores with a ~ 1 and f3 ~ 0.02 should easily fragment. 

More careful analytic work by Tsuribe & Inutsuka (1999a, b; hereafter TI) on 
the dynamics of a collapsing, rotating ellipsoid, including the inward propagation 
of the rarefraction front from the boundary of the cloud, revised the classic 
criterion for the fragmentation of uniform density, rigidly rotating isothermal 
cores. Their numerical simulations included the effect of small white noise, 
bar-mode, and ring-mode density fluctuations. The TI criterion predicts that 
cores in a state of marginal virial stability shoTild never fragment. They report 
their criterion to be in in good agreement with their numerical results. The TI 
criterion clarifies almost three decades of previous numerical and analytic work. 

In summary, the classic a—/3 criterion, which predicted ubiquitous fragmen- 
tation for slow rotators, is now known to be invalid. Moreover, the TI criterion 
predicts ubiquitous single-star formation for realistic cores near virial equilib- 
rium. Significantly, this means that the problem of binary formation, even in 
idealized uniform cores, has remained unsolved to date. 

However, TI's results are limited in that their choice of perturbation spec- 
trum is not physically motivated; in particular, none of the chosen perturbation 
spectra are consistent with the linewidth-size relation observed in molecular 
cloud cores (Larson, 1981). More important, the amplitudes chosen for the 
density perturbations are unrealistically small: typical cores are observed to 
have turbulent linewidths comparable to the isothermal soundspeed (Jijina et 
al. 1999), and so the perturbations present in real molecular cloud cores will 
typically be nonlinear. 

The outcome of the initial collapse phase hinges upon whether the timescalc 
for the overall collapse of the core is greater or less than the timescale for the 
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Figure 1. An isothermal uniform model with a = 1.0, (3 = 0.1, with 
25% m = 2 density perturbation, shown after a compression of over 
seven orders of magnitude in density. Note the formation of a binary 
system. 

growth of perturbations. When the timescale to grow the perturbation is greater 

than that of the overall collapse, a single fragment is formed; in the opposite 
case , a binary or multiple system may be formed. Nonlinear perturbations play 
a key role in this process, since the timescale for the gravitational instability to 
grow a nonlinear perturbation is always less than that for a linear perturbation. 

To test the sensitivity of the TI fragmentation criterion to the choice of 
perturbation amplitude, we apply an m = 2 density perturbation to a variety of 
uniform clouds. We find that a modest increase increase in the amplitude of the 
applied spectrum can cause models significantly removed from their fragmenta- 
tion transition curve to fragment. In figure 1, we show such a model system with 
a = 1.0, P = 0.1, and an applied amplitude of 25%, which produced a binary. In 
contrast, the same background model run with a slightly smaller perturbation 
amplitude of 10% produced a single object at a comparable compression. 

This result illustrates that nonlinear perturbations, which are typically found 
in molecular cloud cores, play a critical role in fragmenting uniform isothermal 
spheres close to marginal virial stability. We turn our attention in the next 
section to a more realistic choice of background core models and applied spectra. 

3. Turbulent Low-Mass Molecular Cloud Cores: Initial Conditions 

The statistics of prestellar versus embedded stellar molecular cloud cores suggest 
that the prestellar phase is greater than the free-fall time of the core; hence, the 
core should be in an approximate state of equilibrium (Beichman et al. 1986, 
Myers et al. 1995, Lee Sz Myers, 1999). Accordingly, the initial conditions for 
a typical model we simulate consist of self-consistent hydrodynamic equilibria 
of turbulent, marginally stable, isothermal molecular cloud cores. In contrast, 
nearly all calculations that appear in the literature deal with clouds unrealisti- 
cally initially far from equilibrium (e.g. Boss & Bodenheimer, 1979, Boss 1991, 
Burkert & Bodenheimer, 1993), leading to highly supersonic collapse velocities. 

We perturb the cloud using Gaussian spectral perturbations on the veloc- 
ity field with a k~^d^k spectrum consistent with Larson's linewidth-size relation 
(Larson, 1981; Dubinski, Nararyan, &; Phillips 1995). A limitation of our present 
procedure is that the density perturbations form only after the velocity pertur- 
bations have time to evolve; by the time the density perturbations have formed, 
the turbulence has decayed somewhat. 

The Bonnor-Ebert values for mass and radius are rescaled to account for 
turbulent support : 

M^M(l + A^2)3/2^ 

R{l + M'^y/'^, (2) 

where A4 = a/cg and a is the ID nonthermal velocity dispersion. We take M 
as the first dimensionless parameter specifying our model cores. 
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As the cloud collapses, the heating rate due to compression rises, while 
the ability of the gas to cool is impeded by greater opacity due to the higher 
column density of dust. As a result, the collapse evolves from isothermal to 
quasi-adiabatic. To follow this evolution in this current set of calculations, we 
use a barotropic equation of state that transitions from isothermal to polytropic 
at a characteristic density of pcrit, i-O-, 

^(p)=4oP + ^stiffP^^*'«, (3) 

where P is the gas pressure as a function of density p, Ciso is the isothermal 
soundspeed, 7stiff is the adiabatic exponent in the optically thick regime, and 
i^stiff is a parameter fixed by the requirement that the isothermal and the adia- 
batic pieces of the equation of state are equal at a characteristic density Pcrit- 

„l-7stifl 

^stiff = (4) 

(Boss et al. 2000; Klein et al. 1999). We take the density ratio Pcrit/Pedge as our 
second dimensionless parameter. As Boss et al. (2000) point out, this barotropic 
equation of state is very approximate, but it enables us to take a first step in 
understanding what role non-isothermal effects play in the more general case. 

Hence, our barotropic models are completely specified in a two-dimensional 
parameter space: {M., Pcrit / Pedgc) ■ Both of these parameters are constrained by 
observation. Nearby star-forming low-mass cores have typical values of ~ 
0.5 — 1.5 (Jijina et al. 1999). We choose a fiducial value of the external density of 
Pedge = 5 X 10"^'' gm cm~^ and a critical density ratio in the range Pait/Pedge ~ 
10^ — 10^. Following Masunaga &l Inutsuka (1998), the critical density ratio can 
be related to the Planck mean opacity Kp by noting that an optically thick, 
gravitationally bound structure at Pcrit will have a characteristic scale of about 
the local Jeans length Aj. Hence, 

KpPcritAj ~ 1 (5) 

With some rearrangement, and assuming we have a cosmic gas mix, we find 

Pedge Pcrit /Pedgc J-U-tV 

Values of appropriate for cold, dense molecular cloud cores are model- 
dependent on dust grain composition, carbon and silicate inclusions, icy man- 
tles, and dust size distribution in the presence of sticking and destructive effects. 
Ranges Kp ~ 0.02 - 0.2 cm^/g are cited in the literature (Preibisch et al. 1993). 
This implies a range of Pcrit/Pedge ~ 10^ ~ 10^- 

Below we explore the effect of turbulence by fixing Pcrit/Pedge = 10^, and 
vary A4 from ~ 0.5 — 3.0. This choice of the stiffening parameter is slightly 
too stiff in comparison to realistic dust grain models. However, the enhanced 
stiffening allows us to evolve our calculations much longer in accretion clock 
time. A future set of calculations will relax this assumption, and will evolve 
with a relatively soft equation of state for comparison. 
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The principal physical processes that wc are neglecting in this study are 
those associated with magnetic fields. Although ambipolar diffusion diminishes 
the importance of the magnetic field relative to both gas pressure and grav- 
ity as the collapse proceeds (Ciolek & Mouschovias, 1994), the role of the field 
in the dynamical processes underlying fragmentation is almost completely un- 
known. However, in order to fully understand the dynamics of the collapse, it 
is worthwhile to carry out careful unmagnetized as well as magnetized studies; 
our current program of research represents a step in this direction. Wc also 
note that some recent studies have suggested that the magnetic field may be 
less important than previously believed (Padoan &; Nordlund 1998). 

4. Results 

We have conducted a preliminary survey of the parameter space of turbulent 
cores by holding the equation of state fixed (i.e., setting Pcrit/Pcdgc = 10^) while 
varying the turbulent Mach number A4 over the range of typical values seen in 
low-mass cores (median Ai ~ 0.5 — 1.5, with the low and high values correspond- 
ing to Taurus and Orion, respectively), to small cluster- forming clumps {Ai ~ 
several) (Jijina et al. 1999). 

All calculations reported utilize our AMR code, which obeys the local Jeans 
condition by resolving the local Jeans length with at least four cells at all points 
in space and time (Truelove et al. 1997; Truelove et al. 1998). Our initial, 
unperturbed, centrally condensed equilibrium isothermal spheres are chosen to 
be very slightly subcritical, with a density contrast of edge to center of 10, and 
a mass of 99.9% that of the critical Bonnor-Ebert value. As a result, without 
turbulence, the models remain stationary indefinitely. In order to obtain physical 
values from our dimensionless models, we set the density at the edge of the cloud 
at Pcdgc = 5 X 10~^^ g cm~^ as remarked above and the isothermal soundspeed 
at Ciso = 1-8 X 10^ cm s~^ (corresponding to T ~ 10 K and cosmic abundances). 
A fiducial turbulent core with = 1.0 then has a radius of i? = 0.07 pc and a 
mass of M = 2.9Mq. The characteristic timescale for the evolution of the entire 
cloud is the free-fall time evaluated at the edge density, tg = (37r/32G/9cdgc)^^^ = 
3.0 X 10^ yr. These fiducial choices are typical of nearby star-forming regions, 
though a substantial variance exists within in each molecular cloud, and from 
cloud to cloud. 

4.1. M = 0.58 

This model corresponds roughly to the median level of turbulence seen in the 
Taurus molecular cloud, which has the lowest level of turbulence of any of the 
nearby star- forming molecular clouds (Jijina et al. 1999). In the fiducial choice 
of edge density and sound speed specified above, this model has R = 0.06 pc 
and M = 1.58Mq. A single fragment is formed, surrounded by a centrifugally 
supported protostellar disk (see Fig. 2). By the end of the calculation, which 
was run for a time ~ 0.5tg after the time of formation of the fragment, the total 
amount of mass in the fragment is roughly half of the initial core mass. This is 
consistent with a nearly constant mass accretion rate over this time interval. 
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Figure 2. The log of the projected column density along the z-axis 
of an = 0.58, Pcrit/Pcdgc = 10^ run. A single fragment has formed, 
and followed to the point it has accreted about 50% of the total mass 
of the core. 

Figure 3. The log of the projected column density along the z-axis of 
an = 1.73, Pcrit/Pedge = 10^ run. A binary system is seen to emerge 
here, which has accreted about 10% of the total mass of the core. 

4.2. M = 1.7 

In the opposite extreme of turbulence seen in nearby molecular clouds is Orion, 
which has a median value A4 ~ 1.5 (Jijina et al. 1999). We ran a calculation 
with turbulent level of support consistent with this value, = 1.7. Using the 
fiducial edge density and sound speed, we find this model has i? = 0.1 pc and 
M = 8.2Mq. It appears that a binary system has formed, which by the endpoint 
of the calculation has accreted about 10% of the total mass of the core. The 
binary system is seen to emerge along the ends of a filamentary structure in 
the turbulent flow (see Fig. 3). Simple estimates of the mass flux through the 
filamentary structure reveal that the density within the filament is typically one- 
scvcral orders of magnitude greater than the ambient medium density, and hence 
that virtually all mass flux occurs along filaments. Both of the members of the 
binary appears to have formed a protostellar disk. Moreover, the disks appear 
to be misaligned relative to one another, as Smith et al. (Smith et al. 1999) 
suggest is the case for T Tau. We note that at this point in the evolution, the 
binary has a very wide separation, ~ 10^ AU. It is possible that this system 
will evolve into a much closer binary, but at present this is too computationally 
demanding to determine. 

4.3. M = 2.9 

Lastly, extending our models into the range of more turbulent clumps that can 
form clusters of stars, we evolved one model with A4 = 2.9. Using the fiducial 
edge density and sound speed, this model has R = 0.16 pc and M = 29M0. A 
complex network of shocked filamentary structures rapidly evolves, and at least 
six fragments emerge from the flow. However, due to computational constraints, 
we have been able to evolve this model only to the point that 0.1% of the total 
mass of the core has been accreted. The fragments appear to form within the 
filamentary structures, at those locations in the flow where filaments cross or 
bend. We are examining the possibility that the mass accretion flow in the 
filaments is largely responsible for the growth of the protostellar fragments for 
such higher Mach cases, though the duration of this calculation is too short to 
draw firm conclusions. 



5. Conclusions 

Beginning with models that have density contrasts, linewidths, radii, masses, 
and projected linewidth gradients consistent with observations of low-mass molec- 
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ular cloud cores, we have demonstrated that turbulence plays a critical role in 
their fragmentation. Cores that are significantly subsonic do not fragment in 
our models, whereas transonic and supersonic cores typically do. Since we in- 
ject turbulence on a level consistent with what is actually observed, our results 
suggest that turbulent fragmentation is likely to play a critical role in the forma- 
tion of binary and multiple systems in molecular cloud cores. Our preliminary 
calculations suggest that filament formation is ubiquitous in turbulent cores. 
The protostcllar cores appear to form within the filamentary structure, and the 
filaments themselves may be the principal conduit for mass transfer onto the 
protostcllar disk. 

However, several important issues remain to be resolved. First, our models 
have difficulty producing binary systems for low turbulent Mach numbers (i.e. 
M < 0.5), whereas observations in Taurus indicate no diminution in the number 
of binaries in this range of M. Second, all models must be evolved for a global 
freefall timescale in order to conclusively answer questions concerning the out- 
come of fragmentation in the cores. Lastly, many different turbulent realizations 
must be used in order to gather statistics regarding final orbital separation and 
mass ratio, which must also be checked against observation. We plan to address 
all of these issues in a subsequent paper (Fisher et al. 2000). This gauntlet of 
tests provides us an unprecedented opportunity to examine the fragmentation 
paradigm in the context of realistic models. 
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